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Abstract 

A  network-theoretic  framework  has  been  developed  to  model  unsteady  flows.  Network  anal¬ 
ysis  has  been  utilized  to  examine  a  wide  variety  of  large-scale  problems,  including  control  of 
disease  transmission,  robust  operation  of  power  grids/internet,  and  uncovering  brain  functional¬ 
ity.  The  present  investigation  extends  the  network-theoretic  approaches  to  describe  the  complex 
nonlinear  behavior  of  unsteady  fluid  flows  with  particular  focus  on  vortex  interaction  and  modal 
energy  transfer.  Three  canonical  unsteady  fluid  flow  problems  are  examined  with  network  anal¬ 
ysis:  (1)  discrete  vortex  dynamics,  (2)  modal  interaction  during  wake  formation  in  cylinder 
flow,  and  (3)  two-dimensional  isotropic  homogeneous  turbulence.  Graph  theory  is  utilized  to 
extract  important  fluid  flow  interactions,  which  enabled  the  development  of  a  sparsified  vortex 
dynamics  model.  The  network-based  approach  to  examine  the  energy  transfer  between  modal 
fluid  flow  structures  has  allowed  for  intuitive  understanding  of  how  perturbation  or  control  input 
can  alter  the  dynamics  of  the  global  fluid  flow.  At  last,  it  was  revealed  that  two-dimensional 
turbulence  has  an  underlying  scale-free  network  structure.  The  findings  from  the  present  study 
lay  out  a  foundation  to  perform  further  network-based  characterization  and  control  to  modify 
the  collective  behavior  of  unsteady  fluid  flows. 
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1  Introduction 

1.1  Objective  and  Accomplishments 

The  objective  of  this  STIR  research  was  to  develop  a  network-theoretic  framework  to  model  un¬ 
steady  vortex-dominated  flows.  To  achieve  the  goals  of  this  project,  we  have  considered  modeling 
and  analysis  of  three  canonical  fluid  flow  problems;  namely  (1)  discrete  vortex  dynamics,  (2)  un¬ 
steady  laminar  flow,  and  (3)  two-dimensional  turbulence.  Within  the  9-month  performance  period, 
we  have  been  able  to  examine  a  discrete  vortex  dynamics  problem  in  detail  and  develop  a  sparsified 
dynamics  model.  Furthermore,  we  were  able  to  make  headways  in  obtaining  preliminary  results 
in  capturing  energy  transfers  in  unsteady  laminar  flows  and  characterizing  turbulent  flow  using 
network  analysis. 

In  particular,  we  have  been  able  to  accomplish  the  following  items: 

1.  Network  representation  of  discrete  vortex  interactions 

Established  network  representations  of  discrete  vortex  dynamics  in  two-dimensional  flow. 
Identified  key  vortex-to- vortex  interactions  through  of  use  of  graph-theoretic  techniques.  Con¬ 
structed  a  sparsified- dynamics  model  based  on  spectral  graph  theory  to  create  computationally 
inexpensive  representation  of  full  nonlinear  dynamics. 

2.  Examination  of  energy  transfer  in  modal  network  of  fluid  flow 

Performed  reduced-order  modeling  of  two-dimensional  flow  over  a  cylinder  using  proper 
orthogonal  decomposition  mode  based  Galerkin  projection.  Modal  network  representations 
are  constructed  to  track  the  energy  flow  between  the  modes.  Feedback  control  was  employed 
to  force  the  flow  to  the  base  state  and  their  implications  on  network  interactions  of  the  modes 
were  examined. 

3.  Characterization  of  two-dimensional  turbulence 

Analyzed  the  network  structure  of  two-dimensional  isotropic  homogeneous  turbulence. 
It  was  found  that  turbulent  flow  was  an  underlying  scale-free  network  that  characterizes  the 
vortical  interactions. 

1.2  Background 

Network  describes  how  components  are  linked  to  one  another.  We  can  represent  the  components 
by  points  (nodes)  and  the  connections  by  lines  (edges)  through  a  mathematical  abstraction.  The 
structure  comprised  of  these  nodes  and  edges  is  called  a  graph,  which  has  been  studied  in  detail 
in  the  field  of  graph  theory  [1].  Network  analysis  is  concerned  with  the  study  of  graphs  as  well 
as  the  interaction  and  evolution  of  the  variables  of  interest  on  graphs  over  time  and  space  [2]. 
The  framework  developed  in  network  analysis  and  graph  theory  can  provide  insights  into  how  the 
structure  of  a  network  can  influence  the  overall  dynamics  taking  place  on  the  network. 

Network  analysis  and  graph  theory  by  nature  are  very  fundamental  and  generic,  which  enable 
them  to  impact  a  wide  range  of  applications,  including  the  analysis  of  biological  and  social  networks, 
study  of  traffic  flows,  and  design  of  robust  power  grids  [2],  Biologists  and  medical  scientists  use 
network  analysis  to  determine  how  electrical  signals  travel  inside  the  brain  and  how  abnormality  in 
the  brain  network  connections  can  affect  the  normal  functionality  [3,  4].  In  epidemiology,  researchers 
model  the  outbreak  of  diseases  on  the  population  network.  Locations  with  high  concentration  of 
population,  such  as  airports,  stations,  schools,  and  hospitals,  can  be  represented  on  a  network  with 
large  number  of  connections.  Identifying  such  locations  is  especially  critical  when  containment 
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measures  are  designed  to  control  outbreaks  of  HIV  [5] ,  S ARS  [6] ,  and  Influenza  [7,  8] .  Each  of  these 
diseases  has  different  dynamics  and  an  associated  network  structure.  Once  the  high-risk  groups 
and  areas  are  identified,  network  analysis  can  assist  in  designing  and  implementing  prevention  and 
combat  strategies  in  the  swiftest  manner  with  limited  resources  [9,  10].  Network  analysis  can  also 
reveal  how  a  group  of  people  are  socially  connected  to  one  another  and  examine  how  subgroups 
within  a  population  are  interlinked  in  a  complex  manner  [11],  Moreover,  network  analysis  has  been 
utilized  in  electrical  engineering  to  determine  the  voltages  and  currents  associated  with  electrical 
circuits  via  graph  representations  [12].  In  the  aforementioned  applications  of  network  analysis,  the 
connections  between  people  or  elements  are  highlighted. 

The  dynamics  of  unsteady  fluid  flow  is  governed  by  nonlinear  interactions.  This  is  especially  true 
for  turbulent  flows  resulting  in  multi-scale  interaction  problems.  While  there  has  been  tremendous 
advancements  made  in  the  areas  of  numerical  simulations  and  stability  analysis  in  fluid  mechanics 
[13,  14,  15],  there  has  been  limited  work  highlighting  the  nonlinear  interactions  responsible  for  the 
collective  dynamics  of  fluid  structures  or  vortices.  In  the  present  work,  utilize  the  network-theoretic 
framework  to  study  unsteady  fluid  flows.  Such  analysis  emphasizes  how  a  collection  of  vortices  or 
flow  structures  influences  each  other  through  a  causal  point  of  view  on  a  network  structure.  We 
believe  that  the  present  study  can  provide  an  alternative  tool  to  analyze  how  vortices  or  flow 
structures  interact  in  the  flow  held  and  support  the  development  of  interaction-based  models  to 
capture  unsteady  vortex  dynamics  or  turbulence. 

In  what  follows,  we  report  the  work  performed  in  three  broad  areas  concerned  with  modeling 
unsteady  fluid  hows  using  network  analysis.  First,  we  introduce  the  framework  of  network  theory 
in  preparation  to  apply  to  the  analysis  of  unsteady  fluid  hows.  We  then  consider  three  canonical 
fluid  how  problems: 

1.  discrete  vortex  interactions  captured  with  sparsihed-dynamics  model; 

2.  energy  transfer  in  modal  network  of  fluid  how;  and 

3.  network-based  characterization  of  two-dimensional  turbulence. 

The  report  concludes  with  some  remarks. 

1.3  Network-Theoretic  Framework 

A  network  is  defined  by  a  collection  of  vertices  joined  by  edges.  These  edges  form  the  bonds  that 
connect  the  vertices  or  the  various  entities  of  the  system  together.  The  edges  can  be  associated 
with  weights  describing  the  importance  of  the  connections  between  the  vertices.  Such  a  graph  is 
called  as  a  weighted  graph.  If  the  edges  of  the  graph  have  an  associated  directivity  (i.e. ,  one  vertex 
influences  the  other  but  not  vice  versa),  the  graph  is  called  a  directed  graph.  If  a  vertex  influences 
itself,  a  self- loop  edge  (edge  reconnected  to  itself)  can  represent  such  an  effect.  Any  undirected 
graph  Q  can  be  described  by  a  set  of  vertices  V  =  {v\,V2,  ■■■vn},  a  set  of  edges  E,  and  a  set  of 
weights  w  associated  with  the  edges,  i.e.,  Q  =  {V,E,w}  [16,  2],  An  example  of  a  weighted  graph 
is  shown  in  Fig.  1  (left)  with  weights  associated  with  the  edges  displayed.  A  complete  graph  K, n 
with  N  =  5,  shown  in  Fig.  1  (right),  has  all  the  vertices  connected  to  each  other.  In  other  words, 
this  graph  has  N  vertices  with  complete  set  of  possible  N(N  —  l)/2  edges  without  any  self-loops. 
We  later  use  a  weighted  version  of  the  complete  graph  to  describe  the  interactions  amongst  a  set 
of  discrete  point  vortices  in  the  context  of  fluid  dynamics. 

Network  connections  can  be  summarized  by  its  adjacency  matrix  Ag  £  M.NxN ,  which  is  given 
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Example  graph 


Complete  graph  K 5 


Figure  1 :  (left)  An  example  weighted  graph;  and  (right)  a  complete  graph  /C5. 


by 


Wij  if  (i,j)eE 
0  otherwise. 


(1) 


In  an  unweighted  graph,  the  weight  is  set  to  unity  if  ( i,j )  E  E.  The  diagonal  entries  of  the 
adjacency  matrix  relate  to  the  weight  of  the  self-connecting  (loops)  edges.  The  adjacency  matrix 
of  an  undirected  graph  is  symmetric.  The  degree  ki  of  a  vertex  i  represents  the  summation  of 
the  weights  of  the  edges  connected  to  it  given  by  ki  =  YljLi  \Ag\ij-  The  adjacency  matrix  for  an 
example  graph  Ag  and  the  complete  graph  A/cs  shown  in  Fig.  1  are  given  by 
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where  Q  is  a  weighted  graph  and  JC 5  is  an  unweighted  complete  graph. 

Another  important  matrix  in  graph  theory  is  the  graph  Laplacian  matrix  Lg  E  WNxN , 
is  given  by 


[Lg\ij 


(ki  if  (i,j)  E  E  and  (i  =  j) 

<  -Wij  if  (i,j)  E  E  and  (i  /  j ) 


[0  otherwise. 


which 


(2) 


The  graph  Laplacian  matrix  can  also  be  deduced  from  the  adjacency  matrix,  Lg  =  Dg  —  Ag,  where 
Dg  E  M.NxN  is  a  diagonal  matrix  with  elements  equal  to  degrees  of  vertices,  Dg  =  diag([fej]^;1).  For 
undirected  networks,  it  is  symmetric  and  positive  semidefinite  (singular).  This  Laplacian  matrix 
is  a  discrete  analog  of  the  negative  continuous  Laplacian  operator  (—V2)  and  is  also  called  as  the 
discrete  Laplacian.  It  is  naturally  defined  by  its  quadratic  form.  If  we  have  a  linear  system  of  the 
form  Lgx  =  b,  for  the  vector  x  E  ,  the  Laplacian  quadratic  form  of  a  weighted  graph  Q  is  given 

by  [17] 

X  LgX  —  ^  ]  Wij  i^Xi  Xj)  .  (3) 


The  discrete  Laplacian  is  a  smoothness  indicator  of  x  over  the  edges  in  Q .  The  Laplacian  quadratic 
form  becomes  large  as  x  jumps  over  the  edges  of  Q .  The  definitions  of  adjacency  matrix  and  graph 
Laplacian  form  the  building  blocks  of  the  graph-theoretic  framework.  The  Laplacian  quadratic  form 
will  be  used  below  to  introduce  the  notion  of  spectral  similarity  of  graphs  and  graph  sparsification. 
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2  Network  Representation  of  Discrete  Vortex  Interactions 


We  apply  network  analysis  to  study  the  dynamics  of  a  collection  of  point  vortices  [18,  19].  In  the 
present  work,  the  flow  is  assumed  to  be  two-dimensional,  incompressible,  and  inviscid.  The  spatial 
domain  T>  is  taken  to  be  infinitely  large  and  the  position  of  the  vortices  in  the  flow  field  is  denoted 
by  r  =  (x,  y )  G  V.  We  consider  a  collection  of  N  vortices  to  be  in  the  domain  V  resulting  in  the 
vorticity  field  of 

N 

u{r)  =  ~rj),  (4) 

3= 1 

where  Kj  and  r?  are  the  strength  (circulation)  and  position  of  the  j-th  vortex,  respectively,  and 
5(-)  is  the  Dirac  delta  function.  The  motion  of  the  vortices  is  described  by  the  Biot-Savart  law 

drj 

dt 

where  k  is  the  out-of-plane  unit  normal  vector.  To  determine  the  trajectories  of  the  vortices,  we 
numerically  integrate  the  above  equation.  This  equation  is  the  basis  of  the  discrete  vortex  methods 
that  are  used  to  simulate  unsteady  vortical  flows  in  place  of  discretizing  the  two-dimensional  Euler 
equations. 

The  objective  of  the  present  work  is  to  describe  the  interactions  amongst  a  set  of  point  vortices 
with  graph  theory  and  network  analysis.  The  dynamics  of  the  vortices  will  be  modeled  on  a  sparsi- 
fied  graph  in  place  of  the  full  (complete)  graph  to  derive  a  sparsified  model  that  accurately  captures 
the  motion  of  the  vortices  while  conserving  physical  properties,  such  as  circulation,  Hamiltonian, 
as  well  as  linear  and  angular  impulse. 

2.1  Representing  Vortex  interactions  on  a  Graph 

A  vortex  is  advected  by  the  induced  velocity  from  the  other  vortices  and  not  by  its  own  velocity 
as  seen  in  equation  (5).  Each  vortex  is  influenced  by  all  other  vortices,  which  means  that  each  and 
every  vortex  has  a  connection  to  all  other  vortices  without  any  self-loops.  Thus,  the  discrete  system 
of  N  point  vortices  can  be  represented  by  a  weighted  complete  graph,  /Cat.  The  vertices  or  nodes 
of  the  graph  represent  the  discrete  point  vortices  and  the  edge  weights  represent  the  strength  of 
the  connections  between  them. 

The  motion  of  point  vortices  is  influenced  by  the  strengths  (circulation)  of  the  individual  point 
vortices  and  the  distances  between  them.  The  shorter  the  distance  between  the  point  vortices, 
the  higher  is  the  influence  of  the  vortices  on  each  other.  Also,  a  vortex  with  higher  strength  has 
more  ability  to  influence  the  surrounding  vortices  compared  to  that  with  lower  strength.  We  assign 
the  weight  for  the  interaction  between  two  vortices  to  be  dependent  on  the  characteristic  vortex 
induced  velocity  based  on  the  strengths  of  the  two  vortices  and  the  distance  between  them.  Thus, 
this  implies  that  the  edge  weight  for  the  vortex  network  should  be  proportional  to  n/r,  where  k 
is  the  characteristic  strength  of  the  two  vortices  and  r  is  the  distance  between  them  within  the 
network. 

The  distance  between  the  vortices  is  dependent  on  the  location  of  the  two  point  vortices  |r*j  —  r*j| 
and  the  characteristic  strength  can  be  established  by  considering  their  geometric  mean  \ Kt k3  \ . 
This  choice  makes  the  adjacency  matrix  symmetric  which  allows  the  use  of  the  spectral  sparsification 
algorithm  discussed  in  §2.2.  Therefore,  the  adjacency  matrix  for  the  graph  representation  of  a  set 


EKj  k  X  ( Vj  -  Vj) 
2ir  |  ri  —  rA 2 

3= 1 
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of  point  vortices  can  be  defined  as, 


if  (*>  j)  G  E  and  1  ±  3 

0  otherwise. 


(6) 


Note  that  the  diagonal  elements  of  the  adjacency  matrix  are  zero  as  the  self-induced  velocity 
components  are  zero.  This  translates  to  the  graph  not  having  any  self-loops.  The  interaction 
between  the  system  of  discrete  point  vortices  represented  on  a  complete  graph  can  be  sparsified  for 
deriving  the  sparse  vortex  interaction  model. 

In  addition  to  considering  the  geometric  mean  of  the  circulation  for  the  weights,  algebraic 
mean  of  the  circulation  was  also  considered.  The  computation  using  the  algebraic  mean  (i.e., 
2  ( | Kt  |  +  | Kjj ) / 1 r —  Tj  |)  for  the  weights  perform  similar  to  that  with  the  above  geometric  mean. 
While  we  show  that  geometric  mean  sufficiently  captures  the  vortex  interaction,  there  may  be  some 
room  for  optimizing  the  choice  of  weights.  In  what  follows  we  utilize  the  geometric  mean  in  equation 
(25)  for  the  network  weight. 

2.2  Spectral  Sparsification 

The  sparse  approximation  of  dense  graphs  can  provide  numerous  benefits  for  numerical  compu¬ 
tation.  Graph  sparsification  is  useful  for  designing  computationally  efficient  algorithms  and  helps 
in  identifying  representative  edges  and  associated  weights  [20].  Graph  similarity,  which  forms  the 
basis  for  graph  sparsification,  can  be  derived  in  a  number  of  ways  depending  on  the  desired  similar¬ 
ity  properties.  Distance  similarity  of  graphs  [21]  can  be  achieved  by  sparse  graphs  called  spanners 
that  keep  the  same  shortest-path  distance  between  each  pair  of  vertices  as  in  the  original  graph 
while  cut  similarity  [22]  can  be  achieved  by  maintaining  the  weight  of  the  edges  leaving  a  set  of 
vertices  of  the  sparse  graph  to  be  the  same  as  that  of  the  original  graph.  A  much  stronger  notion  of 
similarity  referred  to  as  spectral  similarity  was  introduced  by  [20].  Spectral  similarity  is  closely  tied 
to  the  Laplacian  quadratic  form  of  the  graphs  as  defined  by  equation  (3).  The  concept  of  spectral 
similarity  directly  leads  to  spectral  sparsification  of  graphs;  that  is  to  create  sparse  graphs,  which 
are  spectrally  similar  to  the  original  graph. 

Spectral  sparsification  is  a  more  general  abstraction  than  cut  sparsihers  and  maintains  spec¬ 
tral  similarity  between  the  sparsified  and  original  graphs.  In  particular,  spectral  sparsification  can 
remove  some  of  the  edges  in  the  graph,  while  maintaining  similar  adjacency  and  Laplacian  eigen- 
spectra.  One  of  the  key  features  of  spectral  sparsification  is  that  it  keeps  the  sum  of  the  weights 
leaving  the  vertex  of  a  graph  constant.  A  spectral  sparsifier  is  a  subgraph  of  the  original  graph 
whose  Laplacian  quadratic  form  is  approximately  the  same  as  the  original  graph  [20]. 

Sparsification  involves  the  creation  of  a  sparse  graph  Qs  from  the  original  graph  Q  based  on 
an  approximation  order  of  e.  The  quadratic  form  induced  by  graph  Laplacian  of  Q  is  maintained 
up  to  a  multiplicative  (1  ±  e)  factor  by  spectral  sparsification  [23].  Thus,  the  sparse  graph  Qs  is  a 
(1  ±  e)-spectral  approximation  of  Q.  The  approximation  order  e  can  vary  from  zero  to  unity.  The 
approximation  with  e  =  0  indicates  that  Qs  is  same  as  original  graph  Q  and  none  of  the  edges  are 
sparsified,  while  e  =  1  relaxes  the  quadratic  form  induced  by  the  sparse  graph  to  within  twice  of 
that  induced  by  the  original  graph.  An  approximation  of  e  =  1  leads  to  a  heavily  sparsified  graph. 
Denoting  the  Laplacian  matrices  of  Q  and  Qs  by  Lg  and  Lgs,  respectively,  the  spectrally  sparsified 
Laplacian  satisfies 

(1  —  e)x1  Lgx  <  x 1  Lgsx  <  (1  +  e)xT Lgx  (7) 

at  least  with  probability  1/2  with  large  N  for  all  x  £  [24],  For  the  example  problem  considered 

later,  these  bounds  are  much  tighter.  This  tells  us  that  Lgs  holds  eigenvalues  similar  to  those  of 
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Lg.  These  spectrally  similar  sparse  graphs  are  found  using  the  spectral  sparsification  algorithm 
based  on  sampling  by  effective  resistance  discussed  below. 

Before  we  discuss  how  a  graph  can  be  sparsified,  let  us  first  follow  the  works  of  [1]  and  [25]  to 
establish  an  analogy  of  a  graph  to  an  electrical  circuit.  If  the  entire  graph  is  viewed  as  a  resistive 
circuit,  we  can  define  a  resistance  on  the  individual  edges  e  =  (■ i,j )  of  the  graph.  According  to 
Thomson’s  principle,  the  potentials  and  currents  in  a  resistive  circuit  distribute  themselves  so  as 
to  minimize  the  total  energy  in  the  network  [1],  This  energy  minimization  principle  leads  to  the 
concept  of  effective  resistance. 

Effective  resistance  between  vertices  i  and  j  is  the  potential  difference  induced  between  them 
when  a  unit  current  is  injected  at  one  vertex  and  extracted  at  the  other  [1,  25].  The  effective 
graph  resistance  (also  called  as  resistance  distance)  is  the  sum  of  the  effective  resistance  over  all  the 
pairs  of  vertices  in  the  graph  Q  [26,  27].  Rayleigh’s  monotonicity  law  states  that  pairwise  effective 
resistance  is  a  non-increasing  function  of  the  edge  weights  [28]. 

In  order  to  obtain  an  expression  of  effective  resistance  for  graph  sparsification,  we  orient  the 
edges  of  the  original  weighted  undirected  graph  Q  with  N  vertices  and  M  edges.  We  can  represent 
any  directed  graph  by  a  signed  edge  (e)-vertex  (v)  incidence  matrix  Bg  E  M.MxN  given  by 

1  if  e  E  E  and  v  is  the  head  of  e 

—  1  if  e  E  E  and  v  is  the  tail  of  e  (8) 

0  otherwise. 

The  row  of  Bg  corresponding  to  an  edge  e  =  (i.  j)  is  given  by  (p,  —  q3 ) ,  where  Pi  and  q3  are 

elementary  unit  vectors  in  the  i  and  j  directions,  respectively.  If  the  edge  weights  of  the  graph 

are  represented  in  a  diagonal  matrix  given  by  Cg  E  RMxM,  we  can  express  the  Laplacian  matrix 
based  on  the  incidence  matrix  as 

Lg  =BTgCgBg=  ^  Wij(Pi  -  qj)(Pi  -  qj)T .  (9) 

i,jeE 


The  above  relation  holds  true  only  for  undirected  graphs. 

For  an  edge  corresponding  to  e  =  the  unit  current  is  injected  at  vertex  i  and  extracted 

at  vertex  j.  Thus,  we  set  the  electrical  current  across  the  edge  to  be  (Pi  —  q3 ) .  The  potential 
induced  by  this  current  at  the  vertices  is  given  by  Lg{Pi  —  q;),  where  Lg  is  the  Moore-Penrose 
pseudoinverse  of  the  Laplacian  matrix  Lg  [25].  The  potential  difference  across  edge  e  =  (i,j)  is 
then  given  by  (Pi  —  q3)7  Lg(Pi  —  q3).  Thus,  for  unit  current,  the  effective  resistance  across  edge 
e  =  (L  j)  which  corresponds  to  the  potential  difference  can  be  expressed  as 

[Re\ij  =  (: Pi  -  QjfLQ  ( Pi  -  Qj)-  (10) 


The  above  expression  is  used  for  computing  effective  resistance  of  the  edges  of  the  graph.  The 
sparsification  of  the  original  graph  Q  =  {V,  E,  w}  is  performed  with  an  algorithm  Sparsify  [24]  to 
produce  a  sparse  graph  Q$  =  {V,  E,  w}  where  w  are  the  weights  corresponding  to  the  sparse  graph 
Qs ■  This  algorithm  is  based  on  the  concept  of  effective  resistance  and  yields  a  (1  +  e)  sparse  graph 
Qs ■  This  sparse  graph  contains  a  reduced  number  of  0{N\og(N) / e2)  edges.  The  procedure  for 
sparsification  Qs  =  Sparsify(^)  is  summarized  below. 

We  first  create  a  list  of  the  edges  E  with  the  associated  weights  of  the  original  graph  Q .  The 
adjacency  and  Laplacian  matrices  of  graph  Q  are  constructed  from  equations  (1)  and  (2).  The 


Moore-Penrose  pseudoinverse  of  the  Laplacian  matrix  Lg  is  then  computed.  For  each  edge  in  the 
edge  list,  the  elementary  unit  vectors,  p{  and  qp  and  edge  weights  we  are  obtained.  Next,  the 
effective  resistance  [Re]ij  corresponding  to  each  edge  is  computed  from  equation  (10). 

A  random  edge  e  =  ( i,j )  from  the  edge  list  of  graph  Q  is  chosen  with  probability  pe  proportional 
to  we[Re]ij.  The  edge  e  is  added  to  the  sparse  graph  Qs  with  the  weight  given  by  we  =  we/qpe,  where 
q  =  8iVlog2(iV)/e2.  We  take  integer (q)  number  of  samples  independently  without  replacement  and 
sum  the  weights  if  an  edge  is  chosen  more  than  once.  The  resulting  graph  becomes  the  sparsified 
graph  Qs  that  satisfies  equation  (7)  with  eigenvalues  similar  to  those  of  Q. 

After  the  random  sampling  procedure,  a  sparsified  adjacency  matrix  Ags  is  obtained.  For 
convenience,  we  define  the  ratio  of  the  sparsified  and  original  adjacency  matrix  weights  as  WtJ , 


Wij 


Wij/wij  if  {i,j)€E 
0  otherwise. 


(11) 


This  sparsification  factor  Wij  is  related  to  the  probability  of  an  edge  e  =  (i.  j)  of  graph  Q  being 
sampled.  The  edges  of  the  original  graph  that  are  not  sampled  during  the  random  sampling 
procedure  have  zero  weights  in  the  sparsified  graph.  These  edges  are  cut  during  sparsification.  The 
spectral  sparsification  procedure  also  redistributes  the  weights  of  the  cut  edges  among  the  other 
edges  of  the  sparsified  graph.  Thus,  cutting  of  graph  edges  is  compensated  by  redistribution  of  the 
weights  to  preserve  spectral  properties  of  the  original  graph. 

The  spectral  sparsification  algorithm  described  above  produces  a  (1  ±  e)  expander  graph,  i.e., 
sparsifier  with  strong  connectivity  properties  when  e  >  l/y/N,  where  N  is  the  number  of  ver¬ 
tices  [24],  With  the  concept  of  graphs  and  spectral  sparsification  now  discussed,  we  consider  the 
application  of  network  analysis  to  discrete  vortex  dynamics  in  the  next  section. 


2.3  Sparsification  of  Vortex  Interactions 

In  this  section,  we  construct  sparsified  representation  of  vortex  interactions  by  sparsifying  and 
redistributing  the  weights  of  the  connections  between  the  point  vortices  as  discussed  in  §2.2.  We 
consider  a  configuration  of  N  point  vortices  in  nc  clusters.  This  setup  allows  us  to  examine  the 
effect  of  sparsification  on  a  cluster  of  vortices.  The  edge  weights  represented  in  the  adjacency 
matrix  are  given  by  equation  (25). 

Let  us  first  consider  a  collection  of  N  =  100  point  vortices  with  nc  =  5  with  each  cluster 
containing  20  vortices.  The  clusters  in  this  setup  can  be  clearly  identified  by  the  use  of  graph  clus¬ 
tering  algorithms.  An  algorithm  maximizing  modularity  [29]  can  be  utilized  for  the  identification 
of  clusters.  One  could  also  use  /c-means  for  cluster  identification  similar  to  the  work  performed  by 
[30] .  The  point  vortices  in  a  particular  cluster  are  given  a  normal  distribution  in  space  with  a  radial 
standard  deviation  of  or  =  0.2.  The  strength  of  the  point  vortices,  within  the  individual  clusters 
have  a  normal  distribution  with  the  mean  of  R  =  0.1  and  a  standard  deviation  of  aK  =  0.01.  The 
setup  for  this  configuration  is  shown  in  figure  2.  The  positions  of  the  clusters  are  normalized  by 
the  average  radial  distance  of  the  centroid  of  the  clusters  from  the  geometric  center  of  the  system 
(Ro)- 

The  original  point  vortex  distribution  modeled  as  a  complete  graph  is  shown  in  figure  3  (top 
left).  The  colors  of  the  point  vortices  (network  vertices)  indicate  their  unweighted  degree  (number 
of  edge  connections).  As  the  original  graph  here  is  a  complete  graph  1C n,  each  vortex  is  connected 
to  every  other  vortex.  Thus,  the  unweighted  degree  of  each  vortex  in  this  example  is  N  —  1  =  99. 
The  sparsity  pattern  of  the  adjacency  matrix  for  the  complete  graph  Ajcn  is  shown  in  figure  3  (top 
right).  We  also  note  in  the  figure  the  sparsity  index  S',  defined  as  S  =  nnon.ze]:o/N2  with  nnon-zero 
being  the  number  of  non-zero  weights  in  the  adjacency  matrix  and  N  being  the  number  of  vortices. 
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Cluster  3 


Figure  2:  The  spatial  arrangement  of  N  =  100  point  vortices  with  five  clusters. 


As  the  diagonal  elements  of  the  adjacency  matrix  have  zero  weights,  nnon_zero  =  N(N  —  1)  for  the 
complete  graph.  For  N  =  100,  a  sparsity  index  of  S  =  0.99  is  obtained  for  the  original  graph. 

We  now  perform  spectral  sparsification  on  the  original  graph  (Q  =  /Cat)  using  the  algorithm 
Qs  =  Sparsify(^)  described  in  §2.2  to  obtain  the  spectrally  similar  sparse  representations  of  the 
complete  graph  based  on  approximation  order  e.  The  corresponding  adjacency  matrix,  Ags  can 
also  be  found.  The  vortex  network  for  sparse  graph  and  the  sparsity  pattern  of  Ags  corresponding 
to  approximation  order  of  e  =  0.5  are  shown  in  Fig.  3  (middle).  We  observe  that  the  number  of 
edges  (connections)  between  the  vortices  in  cluster  3  and  5  are  reduced.  Because  the  weights  of  the 
adjacency  matrix  are  inversely  proportional  to  the  distance  between  the  vortices,  a  large  number  of 
edges  between  the  vortices  of  the  clusters  with  larger  distances  are  cut  during  sparsification.  The 
sparsity  index  decreases  from  S  =  0.99  for  the  original  complete  graph  to  S  =  0.741  for  the  sparse 
graph  (e  =  0.5),  reducing  the  number  of  connections  by  approximately  25  percent. 

The  vortex  network  and  sparsity  patterns  for  approximation  order  of  e  =  1  are  shown  in  figure  3 
(bottom).  We  observe  dense  representation  of  the  connections  between  the  clusters  1  and  2.  Thus, 
the  proximity  in  clusters  is  identified  both  in  the  network  structure  and  the  sparsity  patterns. 
Also,  we  realize  from  the  sparsity  patterns  that  the  connections  between  the  vortices  in  a  particular 
cluster  are  maintained  while  majority  of  the  inter-cluster  ties  are  cut  resulting  in  a  sparse  graph. 
The  number  of  connections  in  the  sparse  graph  for  e  =  1  are  reduced  dramatically  by  nearly  60 
percent  with  S  =  0.4006. 

For  weighted  graphs,  the  effective  resistance  is  computed  in  such  as  way  so  as  to  maintain 
spectral  similarity.  The  eigenspectra  (a  and  A)  of  the  adjacency  and  Laplacian  matrices  for  the 
sparse  and  original  graphs  are  compared  in  figure  4.  We  observe  from  figure  4  that  the  spectra 
of  the  sparse  and  original  configurations  are  in  good  agreement  as  expected.  The  second  smallest 
eigenvalue  of  the  Laplacian  matrix  represents  the  spectral  gap  or  algebraic  connectivity  of  the  graph, 
i.e. ,  vertices  connected  by  at  least  a  single  path  [2],  As  observed  from  Fig.  4,  the  spectral  gap  is 
greater  than  zero  for  original  and  sparse  graphs  indicating  that  the  graph  is  connected.  Thus, 
sparsification  preserves  the  connectivity  properties  of  the  original  graph.  The  sparsified  graphs 
slightly  underpredict  the  maximum  eigenvalues  of  the  adjacency  and  Laplacian  matrices  compared 
to  the  original  graph.  The  preservation  of  spectral  properties  is  attractive  if  the  sparsified  network 
is  to  be  used  to  describe  the  dynamics  of  the  vortices.  We  expect  the  overall  motion  of  the  vortices 
to  remain  similar  with  sparsified  representations  of  the  vortex  interactions. 

Sparsification  helps  identify  the  individual  clusters  by  creating  subgraphs  within  a  graph  where 
the  density  of  the  edges  between  vertices  is  much  greater  than  that  outside  it.  The  spectral 
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Figure  3:  Original  and  sparsified  graphs  and  their  respective  sparsity  patterns  of  the  adjacency  matrix  for  N  = 
100  vortices.  The  color  of  the  nodes  of  the  vortex  network  represents  the  unweighted  connections  (number  of 
connected  edges)  of  the  individual  point  vortices.  The  color  in  the  sparsity  patterns  of  the  adjacency  matrix 
indicates  the  adjacency  weights  wi3  for  original  graph  Aq  and  wtj  for  sparse  graphs  ASs.  The  empty  white 
spaces  indicate  sparsity  in  the  adjacency  matrix.  Sparsity  index  is  S'  =  nnon.zem/N2,  where  nn0n-zero  is  the  number 
of  non-zero  weights  (elements)  in  adjacency  matrix  and  N  is  the  number  of  vortices.  The  colored  circles  in  the 
vortex  network  represent  the  individual  cluster  groups. 
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(a) 


(b) 


Figure  4:  Eigenvalue  spectra  of  (a)  adjacency  and  (b)  Laplacian  matrices  of  the  sparsified  graphs  with  e  =  0.5 
and  1  in  comparison  to  the  original  graph. 
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Figure  5:  Sparsity  patterns  of  adjacency  matrix  for  N  =  500  and  2000  for  approximation  order  e  =  1.  The  color 
in  the  sparsity  patterns  of  the  adjacency  matrix  denotes  the  adjacency  weights  wi:  with  the  empty  white  space 
indicating  sparsity. 


sparsification  procedure  leads  to  the  reduction  in  the  number  of  edges  to  0(N  \og(N) /e2)  for  large 
N.  We  examine  the  performance  of  spectral  sparsification  on  the  current  example  for  increasing 
the  total  number  of  vortices.  As  in  the  example  that  was  previously  considered,  we  maintain  a 
constant  number  of  clusters  nc  =  5  but  increase  the  number  of  vortices  in  each  cluster  N/nc.  The 
sparsity  patterns  of  the  adjacency  matrix  for  N  =  500  and  2000  with  approximation  order  e  =  1 
are  shown  in  Fig.  5.  We  observe  that  as  the  number  of  vortices  increases,  the  adjacency  matrix  of 
the  sparsified  graph  becomes  increasingly  sparse.  Similar  to  the  case  with  N  =  100  vortices,  the 
majority  of  the  ties  within  a  cluster  are  maintained  while  a  large  number  of  inter-cluster  ties  are  cut. 
The  sparsity  index  decreases  from  S  =  0.1784  to  0.0707  for  N  =  500  to  2000,  which  is  a  substantial 
amount  of  sparsification.  Let  us  further  show  the  sparsity  index  S  for  increasing  N  with  different 
approximation  order  e  in  Fig.  6.  The  sparsity  index  S  decreases  with  an  increase  in  e  and  the  number 
of  vortices.  The  expected  behavior  of  the  sparsity  index  S  =  0(N  log(N))  /  (N2)  =  0(log(N)/N)  is 
observed  for  larger  N.  The  trend  deviates  from  the  expected  behavior  for  lower  N  as  the  availability 
of  the  connections  for  redistribution  of  the  weights  of  the  sparsified  connections  decreases. 

We  observe  that  the  sparsification  algorithm  provides  us  with  a  heavily  sparsified  model  to 
produce  a  computationally  tractable  representation  of  the  vortex-to- vortex  interaction.  We  further 
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Figure  6:  Variation  of  sparsity  index  S  for  sparsification  with  increasing  number  of  point  vortices.  The  reference 
black  line  represents  log(jV)/jV. 


note  that  the  sparsification  algorithm  can  be  easily  parallelized,  which  is  attractive  to  further  reduce 
the  computational  wait  time.  The  bulk  motion  of  clusters  of  point  vortices  could  be  tracked  effec¬ 
tively  by  tracking  their  individual  centroids.  The  present  approach  may  lead  to  algorithms  similar 
to  fast  particle  summation  methods  where  near  and  far-held  effects  are  taken  into  consideration 

[31]. 

2.4  Sparsified-Dynamics  Model 

The  sparsihed  vortex  interactions  can  be  used  for  analyzing  the  dynamics  of  the  system  of  discrete 
point  vortices.  Spectral  sparsification  can  be  viewed  as  adjusting  the  strength  of  the  vortices  in  the 
Biot-Savart  law  through  the  sparsification  factor  W  given  by  equation  (11).  We  can  incorporate 
the  weights  of  the  sparse  graph  into  the  Biot-Savart  law  as 


drj,  _  ST'  w  Kj  kx(rj-  rj) 

dt  ^  27 r  \ri  —  r7j2 

j=i  1 


(12) 


where  Wtj  is  the  sparsification  factor.  The  above  expression  significantly  reduces  the  amount 
of  computation  compared  to  the  original  Biot-Savart  law,  because  the  sparsification  has  made  a 
significant  number  of  elements  of  W  to  be  zero,  as  seen  in  figures  3  and  5. 

For  evaluating  the  effect  of  sparsification  on  dynamics  of  the  point  vortices,  we  consider  the 
same  example  presented  in  §2.3  with  N  =  100  vortices  as  the  initial  condition.  The  objective  here 
is  to  accurately  capture  the  dynamics  of  the  vortex  clusters  with  the  sparsihed  Biot-Savart  law 
given  by  equation  (12).  The  bulk  motions  of  the  clusters  of  point  vortices  are  well  represented 
by  the  centroid  of  the  individual  clusters.  The  original  and  the  sparsihed  dynamics  are  given  by 
equations  (5)  and  (12),  respectively,  and  are  integrated  in  time  with  the  fourth-order  Runge-Kutta 
method.  The  results  from  the  sparsihed-dynamics  model  are  shown  in  figures  7  and  8. 

For  the  present  analysis,  time  is  non-dimensionalized  by  considering  the  cumulative  circulation 
of  the  vortices  of  the  cluster  (T)  and  average  radial  distance  of  the  centroid  of  the  clusters  from 
the  geometric  center  of  the  overall  system  at  the  initial  time  ( R0 ).  As  the  advective  velocity  of  an 
isolated  cluster  of  point  vortices  of  strength  T  at  Ra  is  given  by  u*  =  T /2ttR0,  the  non-dimensional 
time  can  be  deduced  as  tu* / R0  =  tT/2TTR^.  In  the  current  example,  each  individual  cluster  has  a 
collection  of  vortices  with  their  strength  having  a  mean  of  R  =  0.1  (with  a  normal  distribution). 
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Thus,  for  n  vortices  in  a  cluster,  T  =  nR.  The  error  in  position  of  the  centroid  of  the  clusters  for 
the  sparsified  setting  with  respect  to  the  original  setting  given  by  \rt  —  r\  is  non-dimensionalized  by 
the  average  radial  distance  of  the  centroid  of  the  clusters  from  the  geometric  center  of  the  system 
at  the  initial  time  ( R0 ). 

Let  us  consider  the  case  where  sparsification  is  performed  only  once  to  determine  the  sparsi- 
fication  factor  Wl3  before  initiating  the  time  integration.  The  same  sparsification  factor  is  used 
throughout  the  time  integration.  This  procedure  is  inexpensive  as  sparsification  is  performed  only 
once  a  priori.  The  trace  of  the  positions  of  the  centroid  of  the  individual  clusters  and  their  error 
for  approximation  orders  of  e  =  0.5  and  1  for  a  single  sparsification  are  shown  in  figures  7  (top)  and 
8  (top),  respectively.  Considering  the  number  of  connections  cut  between  the  vortices  for  sparse 
approximations,  the  trajectories  of  the  centroid  of  the  clusters  based  on  sparse  vortex  dynamics 
agrees  reasonably  with  those  from  full  dynamics.  On  comparing  the  errors,  we  find  that  the  errors 
with  approximation  order  of  e  =  0.5  is  less  than  those  with  e  =  1,  which  is  expected.  Despite  the 
Biot-Savart  law  being  nonlinear,  the  present  sparsified  approach  achieves  a  reasonable  agreement 
with  the  full  nonlinear  solution. 

If  the  clusters  of  vortices  are  far  from  each  other,  the  connections  cut  between  their  vortices  do 
not  cause  a  large  change  in  the  dynamics  while  the  connections  cut  between  clusters  close  to  each 
other  affect  the  dynamics  considerably.  As  time  progresses,  the  relative  distance  between  vortices 
change,  resulting  in  an  increase  in  error.  For  tY /2ttR 2  >  1,  we  observe  that  the  errors  in  sparse 
vortex  dynamics  increases  to  ~  0.5i?o.  As  the  adjacency  weights  w^j  and  the  sparsification  factor 
Wij  are  based  on  the  initial  position  of  the  vortices,  the  ratio  does  not  hold  over  large  times.  This 
calls  for  adjustment  of  the  weights  and  periodic  sparsification  to  adapt  to  the  dynamically  changing 
position  of  the  vortices. 

We  thus  consider  resparsification  of  the  graph  representing  the  system  of  point  vortices  pe¬ 
riodically  at  (0.ltT/2irRl)  and  (0.01fT/27ri?o)  for  both  approximation  orders.  The  trace  of  the 
centroids  and  error  for  approximation  orders  of  e  =  0.5  and  1  with  resparsification  are  shown  in 
middle  and  bottom  subfigures,  respectively,  of  figures  7  and  8.  We  observe  that  the  error  decreases 
significantly  with  resparsification.  The  error  in  position  reduces  from  0(1O_1)  to  O(10~2)  with 
resparsification.  We  notice  that  the  centroid  trajectories  based  on  the  sparse  and  original  calcula¬ 
tions  become  increasingly  similar  with  resparsification.  Resparsification  updates  the  sparsification 
factor  periodically  based  on  the  position  and  strength  of  the  vortices  and  decreases  the  error  by  an 
order  of  magnitude.  Thus,  the  nonlinear  evolution  of  discrete  vortex  dynamics  is  well  predicted  by 
the  sparsification  techniques. 

One  of  the  advantages  of  sparsification  is  the  decreased  computational  cost  due  to  increased 
sparsity  of  the  connections  between  the  vortices.  This  could  potentially  lead  to  design  of  faster 
algorithms  based  on  sparsification  strategies.  Let  us  first  evaluate  the  offline  cost  of  computing 
effective  resistance  and  random  sampling  required  for  spectral  sparsification.  The  time  required  for 
computing  effective  resistance  tr  and  time  required  for  random  sampling  ts  is  shown  in  figures  9(a) 
and  (b)  respectively.  All  computations  were  performed  in  MATLAB  on  an  iMac  with  3.4  GHz  Intel 
Core  i7  processor.  We  can  observe  that  for  each  edge,  the  computation  time  of  effective  resistance 
is  0(log(IV)),  hence  requiring  0(N2  log(N))  for  the  overall  effective  resistance  computation.  The 
random  sampling  procedure  takes  0(N2)  time  for  all  the  approximation  orders.  We  compare  the 
computational  time  tj  required  at  each  step  of  numerical  integration  of  the  Biot-Savart  law  for  the 
original  and  sparse  configuration.  We  can  observe  from  figure  9(c)  that  the  original  configuration 
takes  0(N2)  time  equivalent  to  the  number  of  edges  in  the  complete  graph,  while  for  the  sparse 
configuration  the  computation  time  is  reduced  to  0(Nlog(N)). 

There  is  always  a  tradeoff  amongst  the  level  of  sparsification,  the  error  that  appears  in  the 
dynamics,  and  the  associated  computational  cost.  The  sparser  the  network  is,  the  faster  the  com- 
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Figure  7:  Trajectories  of  the  centroids  (left)  and  error  in  position  of  centroids  (right)  of  the  vortex  clusters  based 
on  sparse  graph  (e  =  0.5)  and  original  graph  with  time  for  N  =  100,  R  =  0.1.  (a)  trajectories  and  error  based 
on  sparsification  performed  only  at  the  initial  step  while  (b)  results  based  on  resparsification  performed  at  every 
(O.lir /27rJ?o)  and  (c)  results  based  on  resparsification  performed  at  every  (0.01ir/27r R2a).  The  full  nonlinear 
solution  is  shown  with  dashed  lines. 
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Figure  8:  Trajectories  of  the  centroids  (left)  and  error  in  position  of  centroids  (right)  of  the  vortex  clusters  based 
on  sparse  graph  (e  =  1)  and  original  graph  with  time  for  N  =  100,  a  =  0.1.  (a)  trajectories  and  error  based 
on  sparsification  performed  only  at  the  initial  step  while  (b)  results  based  on  resparsification  performed  at  every 
(O.ltT /2nRl)  and  (c)  results  based  on  resparsification  performed  at  every  (O.OltF /2nRl).  The  full  nonlinear 
solution  is  shown  with  dashed  lines. 
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(a) 


(b) 


(c) 


Figure  9:  The  time  required  for  computation  of  (a)  effective  resistance  and  (b)  random  sampling  for  different 
approximation  orders  for  one-time  sparsification.  The  dashed  lines  indicate  the  expected  trends,  (c)  Time  required 
for  numerical  integration  at  each  step  for  original  and  sparsified  configuration  with  e  =  1. 


putation  can  become  but  may  compromise  accuracy.  To  increase  the  computational  accuracy, 
resparsification  can  be  performed  but  can  introduce  an  additional  computational  load.  We  evalu¬ 
ate  the  computational  cost  for  numerical  integration  until  a  total  time  of  0.1tT/2TrR^.  The  time 
required  for  resparsification  and  numerical  integration  for  the  two  different  resparsification  frequen¬ 
cies  (performed  every  O.ltT /2hRq  and  O.OltT /2itRo)  considered  in  this  work  for  original  and  sparse 
configuration  with  e  =  1  is  shown  in  figure  10.  It  should  be  noted  that  for  the  resparsification 
performed  every  O.ltT /2'kR^)1  resparsification  is  performed  once  for  the  total  time  considered.  For 
resparsification  conducted  every  0.01tT/27ri?Q,  resparsfication  is  performed  ten  times.  We  can  see 
that  the  time  required  for  resparsification  at  every  O.ltT /2i tR°  is  considerably  less  compared  to  the 
original  configuration.  With  increase  in  resparsification  frequency,  though  the  time  required  is  less 
compared  to  the  original  configuration,  there  is  increased  cost  associated  with  resparsification. 

To  illustrate  that  the  trajectory  prediction  using  sparsification  cannot  be  performed  naively,  we 
compare  the  results  from  spectral  sparsification  with  those  based  on  random  edge  removal  (without 
redistribution  of  edge  weights) .  We  keeep  the  number  of  cuts  identical  to  that  achieved  by  spectral 
sparsification  for  fairness.  The  sparsity  patterns  for  random  cuts  with  sparsity  S  =  0.7410  and 
0.4006  and  the  corresponding  Laplacian  eigenspectra  are  shown  in  figure  11.  We  observe  that  the 
spectra  from  the  randomly  sparsified  graphs  grossly  under-predicts  that  of  the  original  graph. 

The  trace  of  the  positions  of  the  centroid  of  the  individual  clusters  and  their  error  for  random 
cuts  with  sparsity  S  =  0.7410  and  0.4006  for  one-time  sparsification  are  shown  in  figure  12.  The 
random-cut  approximation  performs  poorly  as  compared  to  spectral  sparsification  and  the  dynamics 
of  the  centroid  of  the  vortex  clusters  are  not  captured.  As  the  connections  between  the  vortices 
are  randomly  cut,  the  centroids  of  the  clusters  move  slower  compared  to  the  original  configuration 
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Figure  10:  The  time  required  for  numerical  integration  for  a  total  time  of  Q.UV/2-kRI  with  different  resparsification 
frequencies. 


(a) 


(b) 


Figure  11:  Sparsity  patterns  of  the  adjacency  matrix  for  N  =  100  vortices  with  random  cuts  with  sparsity  (a) 
S  =  0.7410  and  (b)  S  =  0.4006.  The  color  in  the  sparsity  patterns  of  the  adjacency  matrix  indicates  the  adjacency 
weights  wi3  for  S  =  0.7410  and  0.4006.  (c)  Eigenvalue  spectra  of  Laplacian  matrices  of  the  randomly  sparsified 
graphs  with  S  =  0.7410  and  0.4006  in  comparison  to  the  original  graph.  The  legend  for  the  sparsity  patterns  is 
similar  to  that  of  sparsity  patterns  in  figure  3. 
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Figure  12:  (a)  Trajectories  of  the  centroids  and  (b)  error  i 
on  random  cuts  with  sparsity  S  =  0.7410  and  0.4006  for . 


position  of  centroids  (right)  of  the  vortex  clusters  based 

=  100,  k  =  0.1. 


as  expected.  This  is  due  to  the  loss  of  induced  velocity  with  the  absence  of  weight  redistribution 
from  random  sparsification.  Spectral  sparsification,  on  the  other  hand,  redistributes  the  weights  to 
prevent  the  loss  of  the  overall  interactions  amongst  the  set  of  point  vortices. 

It  is  well  known  that  fast  multipole  methods  [31]  can  also  reduce  the  0(N 2)  velocity  evaluation 
significantly.  The  particle-box  and  box-box  schemes  reduce  the  computational  complexity  from 
0(N2)  to  0(N  log(N))  and  O(N),  respectively.  As  seen  in  this  section,  spectral  sparsification 
reduces  the  computational  load  to  0(Nlog(N ))  for  computing  the  velocity  of  the  vortices.  Fast 
multipole  methods  approximate  the  effect  of  cluster  of  particles  at  a  certain  distance  by  multipole 
expansions  and  organizing  the  particles  to  a  hierarchy  of  clusters  [32],  In  order  to  compensate  for 
the  sparsified  connections  in  spectral  sparsification,  the  weights  are  redistributed  among  the  other 
edge  connections  to  compute  the  sparsified  dynamics.  In  an  analogy  to  fast  multipole  methods, 
the  interaction  list  depends  on  the  connections  that  have  direct  impact  on  the  vortices  and  compu¬ 
tational  accuracy  may  degrade  when  the  vortices  move  in  space.  Similar  to  reconstruction  of  tree 
data  structures  in  fast  multipole  methods,  performing  resparsification  periodically  increases  the 
accuracy  by  re-evaluating  the  associated  weights  as  the  positions  of  the  vortices  evolve  over  time. 
We  note  that  even  with  sparsification,  invariants  of  discrete  vortices  are  conserved  as  discussed 
below. 

2.5  Conservation  Properties 

For  a  discrete  set  of  point  vortices,  angular  impulse  and  linear  impulse  are  among  the  conserved 
quantities.  The  coordinates  of  center  of  vorticity  obtained  from  linear  impulse  and  the  length  of 
dispersion  about  the  center  of  vorticity  obtained  from  angular  impulse,  when  the  strength  of  vortices 
is  of  the  same  sign,  remain  constant.  The  center  of  vorticity  (X,Y)  and  length  of  dispersion  of 
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vorticity  D  are  given  by 


Y  z2i=  1  Y  £i= 1 

-v  \  '.V  ’  1  N  ’ 

2-/j=i  Kj  zJj=i  Kj 

(13) 

j-,2  Eiii - *)2 + (yt- m 

K- 

2^  =  1  K3 

(14) 

where  Ki  is  the  strength  of  the  z-th  vortex  and  (x,  y )  are  the  vortex  positions  [33,  19].  The  circulation 
for  n  vortices  in  a  cluster  is  T  =  nR.  The  total  circulation  of  the  system  of  nc  =  5  clusters  is  given 
by  T*  =  Y.iLi  Kt  =  ncT.  Another  conserved  quantity  is  the  Hamiltonian.  The  Hamiltonian  is 
defined  as  the  interaction  energy  of  the  system  of  point  vortices  and  is  given  by 

1  N 

H  =  ^2  KiKj  log  \J ( Xi  -  Xj )2  +  (yt  -  yj)2.  (15) 

We  non-dimensionalize  the  Hamiltonian  using  the  total  circulation  of  nc  clusters  and  the  average 
radial  distance  of  centroid  of  the  clusters  at  initial  time  R0. 

We  compare  the  error  in  the  Hamiltonian  (He(r / Ra)  —  H(r/R0))/T 2,  the  center  of  vorticity 
{{Xe  —  X)/R0,(Ye  —  Y)/Ra),  and  the  square  of  length  of  dispersion  of  vorticity  (D2  —  D2)/R2 
for  approximation  orders  of  e  =  0.5  and  1  over  time  in  figure  13.  The  preservation  of  center  of 
vorticity  and  length  of  dispersion  implies  the  conservation  of  linear  and  angular  impulse.  Here, 
the  variables  with  subscript  e  denote  those  based  on  a  sparsified  model.  We  observe  that  the  error 
in  the  Hamiltonian  is  of  0(1O~4)  for  e  =  0.5  and  0(1O~3)  for  e  =  1.  The  error  in  the  square  of 
dispersion  is  of  0(1O~6)  while  the  error  in  the  center  of  vorticity  is  of  0(1O~8)  with  sparsification. 
While  not  shown,  resparsification  of  these  vortices  periodically  performed  at  every  (O.ltT /2ttR2) 
and  (O.OlfT /2irR2)  maintain  similar  error  levels  for  the  invariants.  We  also  note  that  the  circulation 
of  each  vortices  has  not  been  altered,  which  conserves  the  individual  and  overall  circulation  over 
time. 

As  seen  in  §2.2,  sparsification  by  effective  resistance  is  based  on  the  energy- minimization  prin¬ 
ciple.  The  effective  resistance  is  defined  such  that  the  total  energy  of  the  system  remains  constant. 
This  is  reflected  well  in  the  conservation  of  the  Hamiltonian  representing  the  interaction  energy  of 
the  discrete  point  vortices.  In  addition,  sparsification  preserves  other  invariants  including  linear 
and  angular  impulse  of  the  discrete  set  of  point  vortices.  We  note  in  passing  that  even  with  respar¬ 
sification,  the  invariants  were  observed  to  be  conserved.  The  conservation  of  the  invariants  for 
the  sparse  configurations  further  encourages  the  applicability  of  sparsification  strategies  on  discrete 
vortex  dynamics. 

2.6  Sparsified-Dynamics  Model  and  Reduced-Order  Model 

We  comment  on  the  similarities  and  differences  between  the  present  sparsified-dynamics  model 
and  the  reduced-order  model.  Both  models  share  the  general  objective  of  deriving  a  model  that 
captures  the  full-order  physics  in  a  distilled  manner.  Reduced-order  models  achieve  such  goals  by 
reducing  the  dimension  of  the  state  variable.  One  such  technique  is  the  Galerkin  projection  of  the 
Navier-Stokes  equations  using  a  set  of  spatial  bases  [34,  35],  such  as  those  determined  from  POD 
[36,  37].  One  can  further  consider  incorporating  the  effect  of  input  and  output  dynamics  by  using 
balanced  truncation  or  the  eigensystem  realization  algorithm  [38,  39].  The  resulting  reduced-order 
models  are  generally  described  by  ordinary-differential  equations  for  the  temporal  coefficients  with 
reduced  dimensions. 
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Figure  13:  Error  in  the  (a)  Hamiltonian,  (b)  square  of  length  of  dispersion,  (c)  re-position  of  center  of  vorticity  and 
(d)  y-position  of  center  of  vorticity  based  on  single  sparsification  for  N  =  100,  R  =  0.1  and  e  =  {0.5, 1}. 


On  the  other  hand,  the  present  sparsified-dynamics  model  does  not  reduce  the  dimensionality 
of  the  state  variable.  We  instead  focus  on  reducing  the  number  of  connections  used  to  capture 
the  overall  dynamics  of  the  fluid  flow.  The  reduction  in  the  number  of  connections  is  performed 
based  on  the  concepts  from  network  analysis  and  graph  theory.  Although  we  cut  a  large  number 
of  the  edges  in  the  graph  representation  of  the  dynamical  interaction,  we  redistribute  the  weights 
associated  with  the  edges  to  maintain  properties  of  the  graph.  This  procedure  is  based  upon 
approximating  the  spectral  properties  of  the  graph  and  does  not  require  selecting  the  spatial  basis 
functions  unlike  the  Galerkin-projection  based  models.  One  nice  feature  of  the  sparsified-dynamics 
model  is  its  ability  to  conserve  physical  variables  such  as  the  Hamiltonian,  angular  impulse,  and 
linear  impulse,  as  discussed  previously.  Hence,  the  sparsified  dynamics  is  able  to  predict  the  full 
dynamics  as  demonstrated  by  the  example  with  discrete  vortices.  By  highlighting  the  interactions 
amongst  a  set  of  vortices,  we  are  able  to  determine  which  interactions  amongst  the  vortices  are 
important  in  guiding  the  overall  motion  of  vortices. 

We  believe  the  sparsified-dynamics  model  has  promising  potential  to  model  various  types  of 
fluid  flow  by  considering  the  modal  structures  as  a  abstraction  of  graph  nodes.  The  application 
of  sparsified  dynamics  models  towards  flow  control  problems  may  also  become  fruitful  as  the  the 
computational  time  necessary  to  capture  the  complex  behavior  of  the  flow  is  reduced  and  the 
interaction  of  nodes  or  flow  structures  are  well-captured,  which  has  been  a  lacking  feature  in  linear 
dynamics  model  and  linear  stability  analysis.  We  anticipate  that  the  presently  proposed  model  can 
lead  to  potential  feedback  control  of  fluid  flow  [40,  41]  but  with  nonlinear  interactions  emphasized. 
One  of  the  open  questions  in  extending  the  sparsified  dynamics  model  for  a  wide  variety  of  fluid 
flow  problems  is  the  choice  of  the  variables  or  modes  to  be  used  for  graph  nodes.  This  issue  is 
currently  being  examined  and  will  be  reported  in  upcoming  studies. 
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3  Energy  Transfer  in  Modal  Network  of  Fluid  Flow 


In  the  above  discrete  vortex  dynamics  problem,  we  examined  the  collective  behavior  of  a  set  of 
potential  vortices  using  network  analysis.  Next,  let  us  study  the  continuous  representation  of  fluid 
flow  using  modal  analysis  with  Galerkin  projection.  In  this  section,  we  consider  the  example  of  two- 
dimensional  incompressible  flow  over  a  cylinder  at  a  diameter-based  Reynolds  number  of  Re  =  100. 
We  construct  an  empirical  Galerkin  model  based  on  POD  (Karhunen-Loeve)  expansion  [42]  of  the 
snapshots  of  data  gathered  from  direct  numerical  simulation  of  incompressible  flow  over  a  cylinder 
using  the  immersed  boundary  projection  method  [43,  44]  as  shown  in  figure  14.  The  work  presented 
below  is  founded  on  the  Galerkin-based  model  of  Noack  et  al.  [45,  46,  47]. 
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Figure  14:  Unsteady  cylinder  flow  at  Re  =  100  and  its  POD  modes. 


3.1  Modal  Energy  Interaction  Analysis 

To  characterize  a  continuous  representation  of  the  flow  field  on  a  network,  we  take  the  finite  Galerkin 
approximation  of  velocity  field 

N 

u(x,  t )  «  u (7V)  =  uq{x)  +  ^2  ai(t)<t>i(x)  (16) 

2—1 

where  uq  is  the  mean  velocity  field,  (pi  are  the  POD  modes,  and 

Oj  =  (u  -  u0,  (pi)  (17) 

are  the  generalized  Fourier  coefficients.  Here,  (a,  b )  denotes  the  inner  product  on  the  computational 
domain.  Mode  pairs  with  alternating  symmetry  properties  can  be  computed  from  proper  orthogonal 
decomposition  (POD)  and  the  energy  of  the  modes  are  given  by 

A i  =  {{u- u0,<pi)2)  =  {a})  ,  (18) 

where  (•)  represents  time  average.  The  kinetic  energy  of  the  7th  mode  is  given  by  Ki  =  Aj/2  and  the 
variation  in  the  kinetic  energy  is  given  by  K  =  \  (( u u'))  computed  from  the  velocity  fluctuations 
u'  =  U  —  Uq. 

The  Galerkin  projection  of  the  incompressible  Navier-Stokes  equations1 

i  1 

—  +  V  •  (mi)  =  -Vp+  — V2u  (19) 

at  Re 

1  incompressibility  is  implicitly  satisfied  by  the  POD  modes  already. 
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onto  the  POD  modes  result  in  a  reduced-order  model  equation  given  by 


where 


N  N 

di  —  Pi  P  ^  ]  Lij o,j  T  ^  ]  CijkUjCik,  i  —  1,2,...  N , 

1= 1  j,k=  1 

T)  =  V2U0)  -  (0i,  V  ■  («0«o)) 

£*1  =  ~  V  ’  “  (&’ V  '  (^1“°)) 
Cijk  =  ~{4>i )  ^  '  (rf’j&k))- 


(20) 


(21) 


Here,  Pi,  Lij  and  Cijk  are  the  constant  mean  shift  coefficient,  linear  Galerkin  coefficient  term,  and 
quadratic  Galerkin  coefficient  term,  respectively. 

As  we  are  interested  in  the  interactions  between  these  modes,  we  track  the  energy  pathways 
between  modes  within  the  flow.  We  use  the  finite  time  thermodynamics  (FTT)  formulation  [46, 
47]  to  deduce  constitutive  equations  for  energy  and  average  modal  amplitude.  Using  Reynolds 
decomposition,  we  can  split  the  modal  amplitude  as  cq  =  (oj)  +  a),  where  m,  =  (a*)  and  Et  = 
7j  <^(a')2).  The  constitute  equations  are  given  by 


N  N  N 

0  =  Pi  p  22  Lijmj  +  22  C ijkmjmk  p  '22  ^CijjEj 

1= 1  j,k=  1  1=1  ^ 

=  Qi  P  Pi 


There  are  primarily  two  energy  interactions  taking  place  between  the  modes.  One  is  external 
energy  Qi  (interaction  with  the  mean  flow)  and  the  other  is  internal  energy  T%  (transfer)  between 
the  modes.  External  energy  is  the  primary  source  of  energy  in  the  flow  derived  from  the  mean  flow. 
External  energy  is  given  by 


N 


Qi  —  2  Ljj  p  2  '  (Cjjj  T  Ciji)Tfl‘j 


1=1 


Ei 


=  2 


1 

L  Re 


N 


(23) 


(<£*,  V20i)  -  v  •  («o0i))  -  (0i,  v  •  (0i«o))  +  ^2{Ciij  P  Ciji)mj 

l=i 


Ei 


The  external  energy  interactions  consist  of  energy  diffusion  (first  term),  convection  (second  term) 
and  production  (third  term). 

The  internal  triadic  interaction  terms  (Fick’s  law)  are  given  by 


Tijk  -  «£ ijk  VKiKjKk  (l  Ri  +  R‘+  J<f  )  ,  (24) 

where  Qjk  —  c ( I Ejjk |  T | Pj.kj |  T | P jik I T | Pjki. |  T | Ckij |  T I Ckji | )  wi 1 1  i  CJijk  —  (0V,  P  dQjQQpQi,  Pp;jk)) 
and  the  parameter  a  is  required  for  balancing  the  overall  kinetic  energy.  Here,  pressure  fluctuations 
are  decomposed  as  P  =  '^QIjkPjkajak-  The  internal  triadic  interaction  energy  can  be  computed  as 
Tj  =  Tijk-  The  triadic  interaction  terms  are  energy  conserving,  i.e.,  YliLi  T%  =  0- 
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3.2  Network  Representation  of  Modal  Energy  Interaction 


We  can  create  network  representation  of  the  modal  energy  interaction  discussed  above.  The  vertices 
or  nodes  of  the  graph  represent  POD  modes.  The  edges  between  the  nodes  indicate  the  transfer 
of  energy  between  the  modes  with  the  weights  denoting  the  magnitude  of  energy  flow.  Also,  the 
orientation  of  the  edges  denotes  the  direction  of  flow  of  energy. 

As  geometric  symmetry  conditions  lead  to  equipartition  of  energy,  we  combine  pairwise  modes 
in  the  flow,  i.e,  {((/q,  </>2),  (</> 3,  04),  Thus,  external  and  internal  energy  interaction  for  mode 

pairs  1  and  2  are  given  by  Q[i,2]  =  Qi  +  Q-2  and  Tn2]  =  Tj  +  T2,  respectively.  Also,  as  the  triadic 
interaction  terms  are  energy  conserving,  we  can  reduce  the  number  of  edges  on  the  graph.  Based 
on  the  sign  of  the  internal  energy  transfer  term,  we  can  classify  the  mode  pairs  into  donors  and 
recipients.  As  donor  mode  pairs  transfer  energy  to  the  recipient  mode  pairs,  negative  sign  of  the 
transfer  term  indicates  donors.  Thus  the  network  representation  consists  of  mean  flow,  donor  nodes 
and  recipient  nodes. 

Based  on  these  observations,  we  can  construct  the  adjacency  matrix  for  the  modal  graph  rep¬ 
resenting  energy  flow  Qe  with, 


iAGe]ij 


Q\n  ,J2]  if  *  =  0  and  jx  >  i,j2  >  i 

<  ~T[h,h]  ^  i  =  1  and  h  >  hh  >  i 
0  otherwise 


(25) 


where,  j\  and  j2  are  mode  pairs.  In  the  above  notation,  index  0  indicates  the  mean  flow  and  1 
indicates  donor  node. 


3.3  Application  to  Linear  to  Nonlinear  Transitional  Flow 

In  this  section,  we  apply  the  above  discussion  to  study  the  growth  of  perturbation  energy  in  base 
state  for  cylinder  flow.  We  study  the  transition  of  fluid  flow  from  the  base  state  to  nonlinear  limit 
cycle.  As  the  problem  considered  (incompressible  flow  at  Re  =  100)  is  unstable,  we  compute  the 
unstable  steady  state  (base  state)  using  the  selective  frequency  damping  technique  [48]. 

Starting  from  the  unstable  steady  state,  we  collect  snapshots  of  the  velocity  field  at  every 
At  =  0.1525  during  which  time  small  numerical  perturbations  are  monitored  to  grow.  This  interval 
is  chosen  based  on  spectral  analysis  (not  shown)  performed  in  the  nonlinear  regime  to  collect  a 
minimum  of  40  snapshots  per  period.  Small  fluctuations  of  the  coefficient  of  lift  Cl  is  observed  near 
the  unstable  steady  state.  As  time  increases,  these  fluctuations  increase  and  they  saturate  in  the 
nonlinear  regime.  The  growth  of  these  fluctuations  from  the  unstable  steady  state  to  the  nonlinear 
(limit  cycle)  regime  is  shown  in  figure  15(a).  We  collect  snapshots  over  Np  =  16  periods.  Based 
on  previous  literature  [49],  we  choose  number  of  modes  N  =  8  for  our  study.  The  corresponding 
variations  in  the  kinetic  energy  of  the  modes  over  time  are  given  in  figure  15(b). 

Network  representations  of  energy  interactions  between  the  modes  at  t  =  30  (linear),  60  (tran¬ 
sitional),  90  (nonlinear  limit  cycle)  are  shown  in  figure  16.  Directionality  of  energy  transfer  is 
indicated  by  arrows  in  the  graph.  We  can  clearly  see  that  mode  pairs  1  and  2  form  the  donor 
modes  of  the  graph.  As  the  windowed  over  which  POD  modes  are  extracts  shifts  over  time  in 
the  present  analysis,  the  POD  modes  changes  along  with  the  mean  flow  for  the  shown  three  flow 
regimes.  The  legend  indicates  the  magnitude  of  the  interactions  in  log  scale.  Also,  we  can  observe 
that  with  increase  in  time,  the  magnitude  of  interaction  with  the  mean  increases,  particularly,  the 
energy  transferred  from  the  mean  to  the  donor  modes.  Also,  as  the  flow  approaches  the  limit  cycle, 
the  contributions  from  and  to  the  mean  flow  are  nearly  equal.  The  weights  of  the  interactions  are 
shown  adjacent  to  the  edges  of  the  graph. 
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Figure  15:  Variation  of  (a)  coefficient  of  lift  and  (b)  modal  kinetic  energy. 
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Figure  16:  Network  representation  of  modal  energy  interaction  for  cylinder  flow  transition  to  nonlinear  regime 
using  a  windowed  approach  at  (a)  t  =  30,  (b)  t  =  60  and  (c)  t  =  90. 


Using  the  windowed  approach  mentioned  above,  the  POD  modes  and  hence  the  nodes  of  the 
graph  change  in  time.  Also,  the  fixed  point  for  the  above  analysis  is  the  mean  flow  as  opposed 
to  the  steady  (base)  state.  This  makes  the  analysis  and  hence  network  representation  increasingly 
complex.  For  transient  flow,  Noack  et  al.  [50]  introduced  the  normalized  mean  field  correction 
called  the  shift  mode.  This  represents  the  shift  of  short-term  time  averaged  flow  away  from  the 
Karhunen-Loeve  space.  The  shift  mode  is  constructed  as  follows. 

N 

0A  =  («O-U8)-y^(uo-US,0i)0i  ,n„s 

i= 1  (26) 

</>A  =  0a/II0aII> 

where  <f>£  is  the  shift  mode  and  us  is  the  base  state  (unstable  steady  state  for  cylinder  flow  at 
Re  =  100).  Here,  as  the  number  of  modes  are  8,  we  consider  the  shift  mode  as  the  9th  mode  in  the 
governing  (dynamical)  equation. 

We  construct  the  nonlinear  Galerkin  model  including  the  shift  mode  and  integrate  it  in  time. 
We  show  the  phase  portrait  of  variation  in  the  shift  mode  with  increasing  time  in  figure  17(a).  The 
variation  in  energy  of  the  mode  pairs  is  shown  in  figure  17(b)  along  with  total  energy  variation  in 
figure  17(c).  We  can  see  that  the  variation  in  modal  energy  in  cylinder  flow  transition  with  the 
Galerkin  model  follows  a  similar  trend  to  the  windowed  approach  but  the  transition  is  delayed. 
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These  results  agree  well  with  those  obtained  by  Noack  et  al.  [50].  As  the  transition  is  delayed 
compared  to  the  windowed  approach,  we  present  the  network  representation  of  energy  interactions 
between  the  modes  at  corresponding  times  as  figure  16  for  the  Galerkin  model,  i.e.,  t  =  90, 195,  270. 
A  shift  mode  is  added  to  the  graph  as  an  additional  node.  The  POD  modes,  mean  and  shift  mode 
and  hence  the  nodes  of  the  energy  graphs  do  not  change  with  time.  The  modal  interactions  are 
similar  to  the  windowed  approach.  We  can  observe  an  additional  interaction  from  the  shift  mode  to 
the  mean  flow  especially  in  cylinder  transition.  This  interaction  however  reduces  as  the  flow  reaches 
limit  cycle  behavior.  As  the  POD  modes  change  in  the  windowed  approach,  evolution  of  nodes  of 
the  graph  becomes  an  additional  concern.  The  Galerkin  model  approach  reduces  the  complexity 
of  problem  by  fixing  the  nodes  of  the  graph.  Addition  of  the  shift  mode  as  a  node  causes  the  fixed 
point  of  the  model  to  be  the  correct  base  state  as  opposed  to  the  mean  flow  for  the  windowed 
approach. 
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Figure  17:  (a)  Phase  portrait  of  shift  mode  in  cylinder  flow  transition.  Variation  in  (b)  modal  kinetic  energy  and  (c) 
total  kinetic  energy  for  galerkin  model  with  shift  mode. 


Figure  18:  Network  representation  of  modal  energy  interaction  for  cylinder  flow  transition  to  nonlinear  regime 
using  a  shift  mode  at  (a)  t  =  90,  (b)  t  =  195  and  (c)  t  =  270. 


3.4  Feedback  Control  of  Modal  Fluid  Flow  Network 

We  implement  feedback  control  on  the  nonlinear  Galerkin  model  and  examine  its  implications  on 
the  modal  energy  interactions.  First,  let  us  observe  the  variation  of  the  Galerkin  coefficients  for 
mode  1,  3  and  shift  mode  9  with  no  control  added  in  figure  19(a).  We  add  proportional  control  to 
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Figure  19:  Variation  of  galerkin  time  coefficients  of  mode  l(oi),  mode  3 (a3),  and  shift  mode  9(ag)  with  (a)  no 
control,  (b)  proportional  control  added  to  mode  1  and  (c)  optimal  control  (LQR)  pinning  mode  1. 


POD  mode  1  using  proportional  control  based  on  Ziegler-Nichols  tuning  [51]  and  show  the  time 
history  of  temporal  coefficients  in  figure  19(b).  Proportional  control  with  forcing  solely  on  mode 
1  guides  the  flow  to  base  state,  i.e.,  modal  time  coefficients  reduce  to  zero  while  the  shift  mode 
settles  to  a  value  obtained  at  base  state.  We  also  perform  LQR  (linear  quadratic  regulator)  control 
on  mode  1.  For  the  LQR  design,  we  consider  the  linearized  governing  equation  as  the  state-space 
equation  given  by  d?;  =  ^2jLi  ^ijaj  +  where  q  is  the  input.  As  forcing  is  introduced  to  mode  1, 
-Bn  =  1  and  all  other  entries  of  B  are  zero.  Considering  identity  matrices  for  the  LQR  weights,  we 
can  compute  the  feedback  term  q  =  —L*jCij.  The  entries  of  the  row  1  of  L*  are  non-zero  as  control 
is  added  to  mode  1.  We  update  the  linear  term  in  the  governing  equation  as  L  —  L*  for  performing 
LQR  control  on  the  nonlinear  governing  equation,  the  results  of  which  are  shown  in  figure  19(c). 
We  can  see  that  optimal  control  (LQR)  works  very  effectively  and  forces  the  flow  to  base  state. 

Proportional  control  added  to  mode  1  amounts  to  adding  a  force  term  to  the  Galerkin  integrated 
equation  for  mode  1  as  shown  in  figure  20(a).  Network  representation  of  energy  interactions  between 
the  modes  at  three  different  times  t  =  65, 140,  220  during  transition  from  oscillatory  (shedding) 
behavior  to  base  state  are  shown  in  figure  20(b),  (c)  and  (d).  We  can  see  that  the  magnitude 
of  mean  interaction  decreases  as  time  proceeds.  The  force  term  added  to  the  node  of  the  graph 
modifies  its  node  weight  and  hence  the  overall  energy  interaction.  Optimal  control  (LQR),  on 
the  other  hand,  amounts  to  subtracting  an  additional  linear  term  from  the  Galerkin  integrated 
equation.  Network  representation  of  energy  interactions  between  the  modes  immediately  obtained 
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Figure  20:  (a)  Force  term  added  to  mode  1  (proportional  control).  Network  representation  of  modal  energy 
interaction  with  proportional  control  added  to  mode  1  at  (b)  t  =  65,  (c)  t  =  140  and  (c)  t  =  220. 

after  application  of  LQR  control  is  shown  in  figure  21.  We  can  see  that  the  magnitude  of  mean 
interaction  decreases. 

In  this  study,  we  created  (continuous)  network  representation  of  fluid  flow  based  on  the  Galerkin- 
based  model  of  Noack  et  al.  [45,  46,  47].  The  network  representations  aided  in  understanding  of 
energy  transfer  between  the  modes  in  the  flow.  Thus,  modal  analysis  allows  us  to  capture  energy 
transfer  from  one  regime  of  flow  to  another.  Transition  of  flow  over  a  cylinder  from  base  state  to 
limit  cycle  behavior  was  examined  to  understand  network  flows  and  fundamental  physics  involved. 
In  addition,  control  strategies  were  implemented  to  modify  these  energy  interactions  and  force  the 
flow  to  base  state.  We  further  plan  to  use  graph-theoretic  tools  to  extract  information  from  network 
representations  and  modify  the  governing  dynamics.  In  addition,  we  plan  to  extend  the  present 
analyses  to  high- dimensional  fluid  flow  problems. 

4  Two-Dimensional  Turbulence  Network 

The  above  two  network-model  formulations  have  focused  on  discrete  vortex  dynamics  and  unsteady 
laminar  flows.  Here,  we  also  apply  network-based  analysis  to  turbulent  flows.  In  particular,  we 
extend  the  Biot-Savart  based  formulation  discussed  in  Section  2  to  two-dimensional  homogeneous 
isotropic  turbulence  [52,  53].  Here,  we  analyze  two-dimensional  turbulence  by  extracting  the  un¬ 
derlying  network  structure  that  captures  the  dominant  vortex  interactions.  Furthermore,  we  show 
that  two-dimensional  turbulence  has  a  scale-free  network  structure  similar  to  other  networks  seen 
in  nature  [54], 
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Figure  21:  Network  representation  of  modal  energy  interaction  for  optimal  control  (LQR)  after  one  period  of 
oscillation. 


4.1  Network  Analysis  Setup 

Network-based  modeling  of  vortex  interactions  presented  in  Section  2  was  able  to  extract  the 
dominant  vortical  network  for  a  group  of  discrete  vortices  [55].  This  was  performed  by  constructing 
a  model  based  on  the  Biot-Savart  law  by  observing  that  a  vortex  influences  the  motion  of  all 
other  vortices  except  for  itself.  The  sum  of  the  induced  velocities  determines  how  each  individual 
vortex  advects  in  the  flow  field.  The  strength  of  the  interaction  between  two  vortices  in  the  flow 
field  is  quantified  by  the  induced  velocity  ug.  The  induced  velocity  from  vortex  i  on  vortex  j  is 
ug^j  =  7i/(27rrjj),  where  the  circulation  of  the  z-th  vortex  is  7 j  and  the  distance  between  vortices 
i  and  j  is  77,-  =  \xi  —  Xj |.  Based  on  the  induced  velocity  between  two  vortices,  the  adjacency  matrix 
can  be  defined  as  Atj  =  ^ ( 1 7-i |  +  \l'j\)/r'ij-  This  measure  of  interaction  is  defined  to  be  arithmetic 
average  of  the  induced  velocity  (without  the  factor  of  27r).  The  geometric  average  Aij  =  17*7^1/^' 

has  also  been  considered  providing  similar  results.  The  key  network  structure  was  determined  by 
performing  graph  sparsification  on  the  full  vortex  interaction  network  using  the  above  adjacency 
definition  [55]. 

In  the  present  study,  we  consider  the  use  of  network  analysis  to  identify  the  underlying  network 
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Figure  22:  Interaction  of  fluid  elements  in  two-dimensional  turbulence.  The  vortical  interaction  between  elements 
i  and  j  having  vorticity  w,  and  cu,  can  be  quantified  through  the  induced  velocities  and  ugj^n.  The  velocity 
induced  by  element  i  onto  j  is  ue,i^j  =  ujiS/(2ivrij),  where  S  is  the  area  of  the  fluid  element  (7^  =  UiS).  The 
vorticity  is  shown  in  the  background  as  a  contour  plot. 


29 


Figure  23:  The  connectivity  distribution  shown  over  the  physical  domain  for  two-dimensional  turbulent  flow  (panel 
a).  Vortex  cores  having  a  high  degree  of  connectivity  act  as  hubs  in  the  turbulent  vortical  network.  Depiction 
of  scale-free  network  structure  over  the  Eulerian  node  points  is  shown  with  the  darkness  of  the  network  edges 
corresponding  to  the  values  of  the  adjacency  elements  (panel  b).  A  majority  of  the  key  edges  are  seen  within 
the  vortical  structures  and  amongst  large-scale  structures.  The  corresponding  degree  distribution  plot  exhibits 
scale-free  characteristics  with  P  «  Ck~'y±A'y  (panel  c).  The  inset  shows  the  circos  plot  of  the  turbulent  network 
depicting  rotational  symmetry.  Panels  (a)  and  (c)  share  the  same  contour  level. 


structure  and  characteristics  of  two-dimensional  unforced  homogeneous  isotropic  turbulence.  The 
adjacency  matrix  is  constructed  by  evaluating  the  local  vorticity  (circulation)  of  fluid  elements  in  an 
Eulerian  discretization  of  a  bi-periodic  domain  as  illustrated  in  figure  22.  Eulerian  grid  points  are 
taken  to  be  nodes  on  the  network  and  edges  are  characterized  by  the  vortical  interactions  between 
the  nodes.  The  two-dimensional  turbulent  flow  is  simulated  by  numerically  solving  the  vorticity 
transport  equation  using  a  Fourier  pseudo-spectral  method  [56]. 

4.2  Scale-Free  Characteristics  of  Two-Dimensional  Turbulence 

Based  on  the  adjacency  matrix,  we  can  identify  which  fluid  elements  have  high  degree  of  connec¬ 
tions  (interactions).  Vortices  with  higher  strengths,  corresponding  to  large  circulation,  are  generally 
larger  and  they  affect  smaller  vortices  as  illustrated  in  figure  23(a).  Hence  larger  stronger  vortices 
seen  in  red  act  as  network  hubs.  Note  that  these  strong  vortices  induce  velocity  over  long  dis¬ 
tances.  In  contrast,  fluid  elements  corresponding  to  smaller,  weaker  eddies,  shown  in  blue,  only 
have  influence  in  their  local  vicinity.  The  vortical  interaction  network  exhibits  a  vast  network  mostly 
connecting  the  vortical  hubs  with  localized  connections  appearing  for  the  weaker  and  smaller  eddies 
figure  23(b).  From  preliminary  analysis,  the  degree  distribution  of  the  two-dimensional  isotropic 
turbulence  network  is  found  to  be  well  represented  by  a  power-law  distribution  P  ~  A:-7  with 
7  =  3.0  ±  0.25  as  shown  in  figure  23(c).  This  distribution  exhibits  scale- free  characteristics  for  the 
turbulence  network.  It  is  interesting  to  observe  that  the  scale-free  nature  of  the  two-dimensional 
turbulence  also  shares  commonality  with  other  scale-free  networks  observed  in  various  fields  of 
natural  and  social  sciences.  The  current  observation  reveals  that  a  scale- free  network  can  describe 
two-dimensional  turbulence  comprised  of  vortices  of  varying  scales  and  strengths.  This  behavior 
cascades  down  to  smaller  scales.  Visualization  of  the  degree  distribution  with  a  circos  plot  (panel 
c)  clearly  shows  the  universal  scale-free  features  of  the  turbulence  network,  as  evident  from  the 
near  rotational  symmetry. 

Characterizing  a  turbulent  flow  as  a  scale-free  network  enables  us  to  view  turbulent  interactions 
in  an  intuitive  manner.  Network  science  provides  insight  into  how  components  of  a  large  system 
are  interconnected.  It  is  well  known  from  network  analysis  that  scale-free  networks  are  resilient  to 
random  perturbations.  However,  attacks  concentrated  at  network  hubs  can  affect  network  dynamics 
in  a  detrimental  manner,  as  seen  in  the  examples  of  aircraft  flight  and  Internet  network  operations. 
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This  provides  us  with  intuitive  explanation  that  turbulent  flows  will  be  resilient  against  small-scale 
forcing  while  the  global  behavior  can  be  easily  modified  by  large-scale  forcing  targeting  large  vortex 
cores.  This  refreshing  viewpoint  allows  us  to  provide  mathematical  guidance  into  the  possibility  of 
taming  turbulent  flows,  which  will  impact  a  wide  spectrum  of  engineering  and  scientific  fluid  flow 
applications. 

5  Summary 

In  this  study,  unsteady  fluid  flows  were  modeled  to  highlight  the  vortical  and  modal  interactions 
using  network-based  approaches.  Three  canonical  fluid  flow  problems  were  chosen  to  develop  the 
network-based  framework.  First,  we  examined  in  detail  the  discrete  vortex  dynamics  to  capture 
the  collective  behavior  of  a  large  number  of  vortices.  Graph  sparsification  was  utilized  to  extract 
the  important  vortical  interaction  network  and  yielded  the  sparsified  dynamics  model.  This  model 
was  found  to  accurately  capture  the  global  motion  of  vortices  and  with  significant  reduction  in 
computational  cost.  Second,  the  modal-interaction  network  was  developed  based  on  the  work  by 
Noack  et  al.  to  examine  how  perturbations  and  control  inputs  propagate  over  the  modal  network. 
We  anticipate  that  this  type  of  approach  will  facilitate  the  analysis  of  higher- dimensional  complex 
fluid  flows.  At  last,  network-based  characterization  of  two-dimensional  homogenous  turbulence  was 
performed.  Preliminary  analysis  has  revealed  that  two-dimensional  turbulence  possesses  a  scale- free 
structure.  This  finding  reveals  insights  into  the  robustness  of  turbulence  network  and  how  we  may 
possibly  control  or  modify  the  dynamics  of  turbulence.  The  fundamental  findings  and  development 
of  network-based  framework  from  this  study  should  enable  future  work  in  analyzing,  modeling,  and 
controlling  a  variety  of  fluid  flows  by  leveraging  the  knowledge  of  interactions  within  unsteady  fluid 
flows. 
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